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ABSTRACT 


A  finite-difference  approach  to  the  description  of  the  electro¬ 
magnetic  field  structure  within  the  region  bounded  by  the  earth  and 
some  outer  spherical  boundary  is  studied.  Within  this  medium,  a 
current  sheet  and  a  free  charge  sheet  are  assumed  to  exist  in  a  spherical 
geometry. 

The  medium  is  assumed  to  be  two-dimensionally  inhomogeneous 
and  to  consist  of  a  plasma  gas.  The  earth  is  considered  to  be  perfectly 
conducting,  and  the  field  distribution  at  the  outer  boundary  is  caused 
to  assume  the  form  of  a  system  of  spherical  waves  which  are  propagating 
outwardly. 
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I.  INTRODUCTION 


An  investigation  into  the  nature  of  the  electromagnetic  field  of  the 
earth  produced  by  both  current  and  free  charge  sources  located  in  the 
upper  atmosphere  is  the  basis  of  this  problem. 

The  medium  surrounding  the  earth  is  taken  to  be  generally 
inhomogeneous  and  anisotropic.  The  two  spherical  boundaries  which 
enclose  the  cavity  of  interest  are  the  earth  or  the  earth’s  core  for  the 
inner  boundary  and  the  magnetospheric  boundary  for  the  outer  confine¬ 
ment  of  the  fields.  It  may,  however,  be  desirable  for  a  different  boundary 
to  be  chosen  for  the  outer  limits  of  the  problem. 

At  some  altitude  is  postulated  a  current  sheet  and  a  sheet  of  free 
surface  charge.  Both  will  initially  be  assumed  to  be  positioned  at  the 
same  altitude  in  order  that  the  notation  used  may  be  simplified.  The 
fields  within  the  boundaries  are  mathematically  examined. 
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H.  MATHEMATICAL  DEVELOPMENT 


Mathematical  analysis  is  now  performed  on  the  problem.  Because 
of  the  geometry  of  the  system,  spherical  coordinates  are  used. 

A.  The  Working  Differential  Equations 

The  appropriate  place  to  begin  the  problem  is  the  group  of 
Maxwell's  equations  which,  written  in  the  standard  vector  notation,  are 


VxH  =  J+|f 

(1) 

_  z?  3B 

VxE=-8T 

(2) 

< 

• 

WJ 

it 

o 

(3) 

V  •  D  =  p  . 

(4) 

For  the  purposes  of  this  problem,  it  is  assumed  that  the 
9D 

displacement  current,  —  ,  is  small  enough  compared  to  the  free  current, 

ot 

J,  so  as  to  be  negligible.  This  approximation  is  generally  valid  for  the 
media  under  consideration  if  the  frequencies  used  are  small.  However, 
the  approximation  is  only  valid  within  plasma  regions  and  does  not  hold 
with  the  relatively  thin  air  layer.  The  frequencies  of  interest  are  those 
frequencies  of  the  micropulsations  in  the  range  from  about  0.  001  cycles 
per  second  to  a  maximum  of  1  cycle  per  second.  This  assumption  reduces 
equation  (1)  to 


2 


For  a  medium  which  is  linear  (non-magnetic),  the  magnetic 

— ^  ^ 

flux  density  B  is  related  to  the  magnetic  field  intensity  H  by  the  permea¬ 
bility  constant  of  the  material. 


B  =  pH  . 


If  the  permeability  is  not  a  function  of  time,  then  equation  (2)  becomes 


V  x  E  =  -p 


8H 

at  ■ 


The  time  dependence  of  all  field  quantities  is  assumed  to  be 


of  the  form  e  ,  Thus, 


E(r,t)  =  E(r)e1U,t 


H(r,  t)  =  H(7)e1U)t. 


There  is  no  loss  of  generality  in  this  assumption  since  any  other  forms 
of  time  variations  may  be  constructed  from  the  assumed  mode  by  a 


Fourier  series. 

The  curl  equation  (6)  is  now  of  the  form 

V  x  E  =  -  iuupH .  (7 ) 

The  factor  e1UJt  is  understood  to  be  associated  with  each  of  the  field 
quantities  and  is  henceforth  not  written. 

For  media  which  are  anisotropic,  the  conductivity  is  usually 
expressed  as  a  dyadic.  The  free  current  density  is  then  given  by  the 


vector  operation 


j  —  (j  •  r. . 
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In  spherical  coordinates  this  operation  is  expanded  into  the  three  com¬ 
ponents  as 


J  =  a  E+a-E  +a  ,E 
r  rr  r  r0  8  r 0  0 


J  =  on  E  +artr,E.  +  a^,E, 
0  0r  r  00  0  00  0 


J  =  a  E  +  a  ,  dE  +  a  ,.E  , 
< p  0r  r  08  0  00  0 


If  the  curl  operations  are  performed  in  this  coordinate  system,  equations 
(5)  and  (7)  become,  in  component  form. 


1  3 

— 7 — a  f  T7  (H  sin  8)  -  —  —*1  =  a  E  +  g  nE  +  q  E 
r  sin  9  L98  0  90  J  rr  r  r8  8  r 0  0 

1  1  9Hr  9 

7  [sin  8  "  ¥r^r  H0  0  =  a0r  Er+  a08E0+  a00E0 

1  d  dHr 

7  [a7(rH8)  '  “9?]  =  a0rEr  +  a0  8Ee+  a00E0 

1  -8  3E8i 

- 7 — o  To  (E,  sin0)  -  -r—  =  -iaipH 

rsin0L9  8  0  90  J  r  r 

1  9E 

7  [sin  0  ~T$  “  ¥7<rE0M  =  ~ia^Hg 

i  a  3E 

7fc<rE8)-irrr]  =  -i^H0 


(8) 

(9) 

(10) 

(11) 

(12) 

(13) 


These  are  the  general  equations  which  are  made  use  of  in  the  solution 
of  this  problem. 

As  they  stand,  equations  (8)  through  (13)  are  rather  un¬ 
manageable  in  the  three  dimensions.  If  it  is  assumed  that  the  conductivity 
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is  not  functionally  related  to  the  azimuth  angle  <f> ,  a  simplification  may 
be  made  which  allows  a  type  of  separation  of  the  preceding  equations. 
For  this  case  it  is  assumed  that  a  symmetry  exists  which  allows  the 
fields  to  be  separated  into  the  form  of 

E(r,  6,  <£)  =  E(r,  9)eirn<^ 

H(r,  B,<p)  =  H(r,  Q)e"rn(p  . 


When  applied  to  equations  (8)  through  (13),  this  assumption  yields,  with 
the  factor  e  ^  understood- 


1  P  0  *1 

- • — h  To  (H,  sin  9)  -  imHn  =  a  E+oftE.  +  a,E 

r  sm  0  L9  8  $  9j  rrr  r8  0  r$  $ 

-  i(rH  )1  =  E  +  aQQE0+anj.E^ 
r  Lsin  0  r  3r  <j>  j  0r  r  00  0  0<£  <p 

13  3Hr 

7[‘37(rH0)  “  ~3{f]  =  °(P rEr+  a</)0E0+  0<p<pK<t> 

Tibe  [h  (E*  sin  e)  •  imEe]  =  ■  iu^Hr 


(14) 

(15) 

(16) 

(17) 


(18) 

(19) 


The  derivatives  in  equations  (14)  through  (19)  are  now 
expanded,  and  the  terms  are  rearranged  to  give 
a  i  9E 

H  COS  9  E _ - _ 0  ,  im  E 

r  iu)}j.r  sin  0  <j>  iwp.r  38  iuujxr  sin  8  8 


(20) 
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E 

r 


_i_  r  cos  8  1  _ 0 

a  Lr  sin  8  0  r  38 
rr 


im 

r  sin  8 


"  areEe 


ar  4>E<p~\ 


(21) 


9H0  i  9Hr  i 

~dr~  ~  a0rEr  +  a0  8E8  +  a00E0  +  7  ~38~  "  7  H8 


3Ert  .  3E 
_ 6  _  1  _ r  1  . 

3r  r  38  "  r  Ee"  1UJfiH0 


3H 


3r 


±  = 


E  '  ^qoEq“  C7q i  E  -  ir?-A  H  -  —  H 
8r  r  88  8  80  0  r  sin  Q  r  r  0 


3E 

=  icujiH  +  H  -  -  E 

3r  8  r  sm  8  r  r  0 


(22) 

(23) 

(24) 

(25) 


The  radial  fields  are  therefore  expressed  in  terms  of  the  tangential 
fields  and  their  derivatives,  while  the  derivatives  of  the  tangential 
fields  in  the  radial  direction  occur  in  terms  of  the  fields  and  their 
tangential  derivatives  only.  A  method  for  the  calculation  of  the  fields 
now  presents  itself. 

Consider  a  quadrant  of  the  constant  0 -plane,  as  shewn  in 

Figure  1,  partitioned  by  a  grid  of  N  points  in  the  8-direction  by  M 

points  in  the  radial  direction.  Between  level  L  (levels  referring  to  the 

radial  direction)  and  level  L  +  l,  a  boundary  exists  which  carries  a 

surface  current  J  and  a  free  surface  charge  distribution  p  . 

s  s 

If  the  tangential  fields  on  the  inner  boundary  are  completely 
specified,  then  their  derivatives  with  respect  to  the  polar  angle  6  are 


Sj  'S&r* 
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Figure  1 


Quadrant  of  Constant  0 -Plane  ?  -  r  grid 


\  =  0  \  =  1. 0 


Quadrant  of  Constant  0-Plane  r  grid 


8 

specified,  and  equations  (20)  and  (21)  may  be  utilized  for  the  determination 
of  the  two  radial  field  components  at  all  points  on  the  inner  boundary.  Inte¬ 
gration  of  the  equations  (22)  through  (25)  then  yields  the  tangential  fields 
at  the  second  level.  This  process  is  continued  until  the  fields  at  level  h 
are  known.  After  proper  boundary  conditions  are  applied  at  the  discon¬ 
tinuity  between  levels  L.  and  L  +  l,  the  fields  at  level  L  +  l  are  determined, 
and  the  marching  process  is  continued  through  level  M. 

The  method  may,  of  course,  be  employed  for  the  entire 
range  of  polar  angle  0  from  0°  to  180°.  However,  the  summetry  of  the 
earth's  dipole  field  about  the  equatorial  plane  leads  to  a  simplification 
in  that  fewer  points  are  needed  if  the  plane  of  6  =  90°  is  thought  of  as  an 
artificial  boundary.  The  field  quantities  are  then  terminated  at  this 
boundary  as  odd  functions  about  0  =  90°  (in  which  case  the  field  is  zero 
at  the  boundary)  or  even  functions  about  the  equator  (in  which  case  the 
derivatives  of  the  fields  with  respect  to  the  polar  angle  is  zero). 

If  the  radial  component  of  the  electric  field  is  chosen  as 
being  an  even  function  of  9  about  9  =  90°,  then  it  can  be  shown  that 

H  and  E  ,  must  also  be  even  functions  about  the  equator  while  H  ,  E„ 

9  <p  r*  9 

and  H  are  odd  functions.  The  converse  is  also  true.  It  is  also  true 
9 

that  the  integer  m  specifies  a  symmetry  of  the  fields  about  the  pole. 

Some  difficulty  may  arise  in  the  use  of  the  differential 
equations  (20)  through  (25)  at  points  near  the  equatorial  axis  if  equal 


swasaatt,  -if  r. 


raaaiMs'tt.  W 
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increments  of  the  angle  theta  are  used.  In  order  that  more  points  be 
located  near  this  axis,  giving  greater  resolution  in  this  area  of  trouble, 
the  independent  variable  theta  is  changed  to  cos  0.  The  two-dimensional 
space  under  consideration  may  now  be  represented  as  in  Figure  2.-. 

The  transformation  from  angle  increments  to  increments 
based  on  the  cosine  of  the  angle  is  effected  by  the  substitution  of  X  =  cos  6 
into  equations  (20)  through  (25).  Then, 


sin  0  =  v  1-X  , 


2  a 


d\  * 


The  working  equations  now  become 


H  =  - 
r 


Hi)  pr 


%  rr T2  dE. 

± _  E  ,  _ f. 

JTTtf  *  i»tir  <>*■ 


_  E 

e 


,  _  1  p  X  H  1  -  X.^  im 

'r '  °rr  L'rs/Tl2  *'  r  3X'rJT^ 


^—7  H  -  a  nEn-  a  .  E  .1  (27) 

1  _  ^2  9  r  0  0  r(p  $ J 


!f!e  =  p  F  ,  p  8A 

9r  °<pr  r  °tpB  Q  °(p<p  <p  r  d\  r  0 


U2  3Er  1  _  .  H 

~  "5x~  "  i  Ee~ 


(b  im  1 

r  =  '  a0rEr  "  a00E0  '  a00E<*>  "  irJ^z  Hv~rU<P 
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9E 

if  -  iujfiHe+‘ 


—J^=  E  -  -  E 

:JTl 2  r  r  <*> 


Now  equal  increments  of  X.  are  used  instead  of  the  0  increments. 


B.  Finite  Difference  Equations 

The  most  straightforward  manner  for  the  solution  of 
differential  equations  such  as  equations  (26)  through  (31)  on  a  high-speed 
digital  computer  is  the  finite  difference  approach.  The  finite -differencing 
of  equations  (26)  through  (31)  follows  with  the  aid  of  standard  formulae 
found  in  Appendix  A.  The  equations  are  then  written  as 


H  <p,q)  - - -==  E  <p,q)+ AdL 

r  iuniWTTz  *  l0»*r 


(P.  q+1)-  E^(p,q-1) 
2A\ 


ujprvl  - 


^2  H  (P,q+1)-H  (P.q-1) 


-y==|  H0(p,  q)  -  or0E0(p,q)-  a^E^p.q)}  (33) 


H0(p+1,  q)  -  HQ(p-l,q) 


=  VEr(p’q)  +  °<t>  eEe(p,q)+  Wp’q) 

{ 

H  (P#q+1)-H  (p,q-l) 

-^r-  [~ - zsr - ]'rHe(p,q) 


EQ(p+i»q)  -  e0(p-i.  q) 


li.yf  Er(p»  q+i)- Er(p»  q-1) 


-  E0(p,  q)  -  itupH^{p,q) 


H^(p+l,q)-  H^(p-1,  q) 


=  -oerEr(p.q)  -  peeEetp-<5)-vE^<p-'J> 


r-feH^q,^VP-q) 


E  (p+i.ql-E  (p-i,q)  lm 

- - 2A? -  =  iUJti  He<P’  q)  +  'Jfr==£  Er(p>  q) 


-7  E0(p,q)  . 


In  equations  (32)  through  (37)  all  conductivities  are  those  at  the  point 
(p,  q),  the  \  is  and  r  is  r  . 

The  final  form  of  the  difference  equations  follows  from  the 
solution  of  equations  (34)  through  (37)  for  the  tangential  fields  at  the  layer 
p  +  1.  The  results  are  the  following  equations: 


's/I  -  \  „  ,  ..  m 

2iu)pr  A\  *VP’  q"  )  +  f— J  E0(P’  q) 

(i-urv  1  -  X 


(38) 
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Er(p,q)  = 


HQ(p+l,q)  = 


E0(p+1»q) 


H  (P+1,  q) 


E^p+i.q) 


a  Wl-X2 


HJp.q)  - 


rr 


2  a  r  AX 
rr 


H0(p,  q+1) 


HA(p,q-l)+  - £=  EJp.q) 


2a  r  AX 
rr 


U)pr 


n/i-x2 


(39) 


H0(p-1,  q)  +  2Ar  a£rEr(P»  q)+  2Ar  a^E^p,  q) 


2Ar  %  A(P’  q>  '  ArrA\X  Hr(P’  q+1) 


Ar\ll  -X2  „  .  ..  2Ar  . 

tax  r  *Wp*q) 


(40) 


Ar 


Jl  *X2 


Ar-s/l-X2 


*9^’  q> '  Er<p’  q+1>  +  E>  q-1} 


^  E0(p,q)-  2iu¥ArHj#>(p,q) 


(41) 


=  H^(p-1,  q)  -  2aQr  Ar  Er(p,  q)  -  2Ar  aQ8Ee(P»  q) 


_  .  „  .  .  2  im  Ar  __  ,  . 

ZAiawVv-q)-^rH r(P’q) 


(42) 


2Ar  . 

—  H^{p,q) 


E  (p-1,  q )  +  2iu)pAr  H  (p,  q )  +  Zl~^L  E  (p,  q )  EJp,  q  ). 
*  6  rsl  1-X2  r 


_  A' 

r  0 


(43) 
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The  functional  form  of  the  difference  equations  is 

Hr(p.q)  =  f[E^(piq)#  E  (p.  q+1)  ,  E^(p.q-l),  E0(pfq)]  (44) 

Er(p,  q)  =  f[H^(p,  q>  ,  H^(p,q+1),  H^(p,q-1),  HQ(p,  q)  , 

E0(p,q),  E  (p,q)]  (45) 

HQ(p+l,q)  =  f  ^HQ(p-l,q)  ,  Er(p,q),  E0(p,q),  E^(p,q) 

Hr(p,  q+1)  ,  Hr(p,  q-1)  ,  HQ(p,  q)]  (46) 

EQ(p+l,q)  =  f^E0(p>l,q)  ,  Er(p,q+1)  ,  Er(p,  q-1), 

E0(p,q)  ,  H  (p,q)  ]  (47) 

H^fp+i.q)  =  f^H^(p-l,q)  ,  Er(p,  q)  ,  E0(p,q)  ,  E^(p,q) 

Hr(p,  q)  ,  H^(p,q)J  (48) 

^(p+M)  =  f[E^(p~i»q) »  He(p»q)»  Er(p>q)>  E^p.q)]-  (49) 


Since  the  forward  difference  is  used  in  the  initial  step, 
H0(p-1,  q),  E0(p-1,  q),  H^(p-1,  q)  and  E^(p-1,  q)  in  equations  (44)  through 
(49)  are  replaced  with  HQ(p,q),  E0(p,q)  ,  E0(p,q),  H^(p,q)  andE^(p,q), 
respectively,  although  the  exact  equations  change  by  a  factor  of  two. 
Because  of  the  symmetry  conditions  imposed,  the  central  difference 


formula  will  also  be  used  at  X.  =  0  and  \  =  1. 


3 
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Initially  all  tangential  fields  at  the  inner  boundary  are  set 

equal  to  zero  except  H  (1, 1).  From  equation  (38)  it  is  seen  that  H  (1,  q) 

0  r 

is  zero  for  all  q.  Equation  (39)  is  used  to  determine  E  (1,  q)  as  a 

function  of  H  (1,  1  ),  all  other  tangential  fields  on  that  level  being  zero. 

Equations  (40)  through  (43)  are  then  applied  to  yield  the  tangential  fields 

at  all  points  on  the  second  radial  level  in  terms  of  HQ(1, 1),  since  all 

quantities  on  the  right-hand  side  of  these  four  equations  are  either  zero 

or  are  in  terms  of  H  (1, 1). 

0 

With  p  =  2  for  the  second  level,  equations  (44)  and  (45)  yield 

the  radial  fields  E  (2,  q)  and  H  (2,  q)  as  a  linear  function  of  H  (1,  1)  and, 

r  r  fc) 

therefore,  all  fields  at  the  second  level  are  specified  as  a  function  of 
H  (1, 1).  The  process  of  obtaining  the  six  field  components  at  each  level 
is  continued  until  all  fields  directly  below  the  current  sheet  and  charge 
distribution  are  known  in  terms  of  Hg(l,  1). 

In  order  that  as  much  accuracy  as  possible  be  retained  in 
the  calculations,  the  level  below  and  the  level  above  the  discontinuity 
should  be  relatively  close  to  this  level  as  in  Figures  1  and  2. 

A  general  discontinuity  is  represented  in  Figure  3.  Surface 
current  and  surface  charge  are  present  on  the  boundary. 


* 
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Figure  3. 

A  General  Discontinuity 

The  boundary  conditions  for  the  fields  around  such  a  dis¬ 
continuity  are 


n  •  (Hj  -  H2)  =  0 

(50) 

n-  (Ej-I2)  =  Ps/eo 

(51) 

“X(HrH2)  =  T 

(52) 

n  X  (E  -  E2)  =  0  . 

(53) 

The  unit  vector  n  for  this  problem  is  the  radial  unit  vector  a^.  Expansion 
of  these  four  equations  in  spherical  coordinates  gives 


i 

X 

i 
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<t> 2 

or,  in  terms  applicable  to  Figure  2, 

Hr(L+l,  q)  =  Hr(L,  q) 

He(L+I,  q)  =  H0(L,  q)  +  J^fq) 
H#(M.,)  =  H  (L.,)-  Js9(q) 
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>n 


Er<L+l,q)  =  Er(L,q}+  pg(q)/eo 


Ee(L+l,q)  -  EQ(L,q) 

E  (L+l,q)  =  E  (L,q)  , 

for  all  q,  1  <  q  <  N. 

These  expressions  relate  the  six  field  components  at  level 
L  +  l  to  the  six  components  at  level  L.  Since  those  at  level  L.  are  already 
known  in  terms  of  H  (1, 1),  both  tangential  and  radial  fields  are  now 

u 

related  to  HQ(1, 1)  and  the  source  terms  directly  below  each  point  on  the 

U 

discontinuity. 

Symbolically, 

Er(L  +  l,q)  =  f  [hb(1,1)  ,  ps(q)] 

E0(L  +  l.q)  =  £[He(l,l)] 

E  (L  +  l,q)  =  f[He(l,l)] 

Hr(L  +  l,q>  =  £[H0(1. 1)  ] 

He(L  +  l,q)  =  £[H9(l,l),  JS(#,<q>] 

H  (L+l,q)  =  £[h0(1,1),  Js0fq)  ]  . 


The  functions  are  linear  with  respect  to  Hg(l,  1),  pg(9)» 


J30(q)- and 


J  Q(q).  Now,  from  equation  (46)  through  (49),  the  tangential  fields 


at  level  L  +  2  are,  in  functional  form, 


HQ(L  +  2,q)  =  f[H0(L,q)  ,  Er(L+l,q),  E0(L  +  l,q)f 

E.(L  +  l,q)f  H  (L  +  l,q  +  l),  H  (L  +  l,  q  -  1), 
<p  r  r 


H0{L  +  l.q)]  (54) 

EQ(L  +  2,q)  =  f[EQ{L,q),  E^L  +  1,  q  +  1),  Er(L  +  l,  q  -  1), 

Eg(L  + 1,  q)»  H  (L  +  l,q)]  (55) 

H  (L  +  2,q)  =  f[H  (L,q),  Er(L+l,q),  EQ(L+l,q), 

E  (L  +  l,q),  Hr(L+l,q),  H  (L  +  l,q)]  (56) 

E  (L  +  2,q)  =  f[E  (L,q),  H0(L  +  l,q),  Er(L  +  l.q), 

E  (L  +  l,q)].  (57) 


Since 

Er(L+l,q  +  l)  =  f[He(l,l),  Ps(q  +  1)] 

Er(L+l,  q-l)=  f [Hq(1,  1),  Ps(q-1)], 

H0(L  +  2,q)  =  f[H0(l,l),  J^lq),  Ps(q)] 

E0(L  +  2,  q)  =  f  [h0(1,  1) ,  ps(q  +  1),  ps(q  -  1)  ,  Jg  0(q)  ] 
H0(L+2,q)  =f[H9(l,l),  Js0(q).  f>g(q)] 

E^L  +  2,q)  =  f[He(l,l),  Jg^ (q) »  P,(q)]. 


so  that  by  equations  (44)  and  (45), 


IS 


Hr(L+2,q)  =  f[He(l,l),  J  (q),  Ps(q),  J  (q  +  1). 

PS(q  +  1)'  •%(q‘1)'  ps(q  '  1)1  Js9(q)] 

Er(L  +  2,q)  =  f[H8(l.l),  Js9(q),  ps(q).  Js0(q  +  1). 

Ps(q  +  1).  JBe(q-i).  Ps(q-1).  Ja0(q)]  • 

It  can  be  seen  that,  as  more  and  more  steps  are  taken  away 
from  the  sources,  the  fields  at  angular  location  of  X.  become  functions 
of  more  and  more  of  the  source  components,  so  that,  at  level  M  (outer 
boundary)  the  fields  may  be  functions  of  all  the  source  terms.  Thus 

H9(M,q)  =  f[H9(l,l>,  Js0(l).  Jse(2).-.  J,e<N).  J^d). 

V2'-"-’  V”1’  ps(1)'  ’  P.W]- 

and  likewise  for  the  remaining  three  tangential  components.  If  all 
tangential  fields  at  the  inner  boundary  are  zero  except  HQ(1,  2),  then 
similar  results  are  obtained  as 

He(M,q)  =  f[HeU.2).  Js^<n),  Jse(n>.  pg(n)], 
n  =  1,  2,  •  •  •  ,  N. 

In  this  manner  the  outer  tangential  fields  are  found  as 
functions  of  each  of  the  tangential  magnetic  fields  at  the  inner  boundary 
and  the  source  terms  at  each  point  of  the  discontinuity.  The  results 


may  be  noted  as 
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E0(M,  q)  =  AqlHQ{l,l)  +  •••  + 


ANHe(l.N)  +  A^(U)+...+ 


W’  N>  +  V^(1)+--  '  +  aqNJs«<N>+  °qlJS0(1) 

+  •  ‘  ‘  +  a'qNJ^(N)  +  gqlPS(1)  +  ■  '  •  +  eqNPs(N>  <58) 


%  <M'  ^ =  BqlH9(1-  11  +  '  *  '  +  BqNH8(1’  N>  +  BqlH^  (1>  11 +  '  '  ‘  + 

B;NH^(1>N)+  ^qlJs^1)+---+  ^qNJ^(N)  + Vse(1) 

+  •••  +  VqNJs0(N)+  6qlPs(1)+"'  +  6qNPs(N) 


H0(M,q)  =  CqlH0(l,l)+-..  +  CqNHe(l,  Nl+C^H^l,  !)+•••  + 

c'qnV1>N)  +  W1’*'”  +  WN)+Vse(1) 

+  '  •  •  +  SqNJse(N)  +  Vs(1)+  '  ' '  +  n,NPs(N)  (60) 


V M'  <*>  =  DqlHe(1,  *>  +  '  •  •  +  DqNH0<1’  N)+  DqlV‘’  ‘)  +  '  '  '  + 


D»  TH  fi,  N)+  C  ^(1)+*  •  •  +  C  j.(N)+  C'J  Jl) 
qN  <f> '  ®ql  s ’  bqN  s <£'  ’  ^ql  s0x 


+  +  CqNJs0<N)+  VqlPs(1)+'**  +  VqNPs(N)*  (6l) 


In  order  to  conveniently  determine  the  coefficients  in  equa¬ 


tions  (58)  through  (61),  we  shall  have  to  treat  the  source  terms  J  ,  J  _ 

s  o)  s  y 
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and  p  as  unknowns.  Then  each  coefficient  is  determined  by  setting  all 
s 

variables  to  zero  except  the  variable  associated  with  the  coefficient 

desired,  which  is  given  the  value  of  unity. 

For  instance,  to  determine  the  value  of  the  4N  coefficients 

A  B  .,  C  .,  D  ,  ;  q  =  1,  2,  .  .  . ,  N;  all  fields  HA  and  H,  at  the  inner 

ql  ql  ql  ql  89 

boundary  points  are  set  to  zero  except  H.(l,  1)  which  is  unity,  and  all 
sources  are  assumed  zero.  Then 


Aql=Ee(M,q) 

Bqi  =  VM,q> 

Cql-He(M,q) 

Dqi  =  VM'q) 

so  that  the  coefficients  are  determined  from  the  outer  boundary  fields 
produced  by  only  a  unity  Hq(1,  1).  Similarly,  a.^,  y^,  and  Q  ^  are 
found  by  letting  all  tangential  magnetic  fields  at  the  inner  boundary  be  zero, 
and  all  source  terms  be  zero  except  Js^(l). 

In  this  manner  all  the  coefficients  are  determined,  and  the 
outer  tangential  fields  are  found  as  known  functions  of  the  tangential 
magnetic  fields  at  the  inner  boundary  and  the  source  distribution. 

Equations  (58)  through  (61)  may  be  rewritten  in  more 
compact  form  as 


I 


H 


E."‘’):iA,n¥"|U;VU) 

n=l 

+  Vs**”*  +  °;nJs9(n)  +  VsH 
%(M,q)  =  I[BqnHe(1'n)  +  B;nH  (1-n) 

n=l 

+  VJ^(n)  +  YqnJsS<n)  +  VPs<n)] 


El 


(62) 


(63) 


Hi«q)  =  ![c,.H9|1'n,(C;„H*(!' 


n= 


nj 


qn  0 

+  §  J  (n)  +  §'  J  (n)  +  ri  P  (n)l 
qn  s <£'  ^qn  s0x  qnrsx  'J 


(64) 


H^M,q)  =  f  [DqnHe(l,n)+D'nH<Ji(l.n) 

n=l 

+  C  J  ,  (n)  +  Q1  J  (n)  +  V  p  (n)1 
^qn  s <p  ^qn  s0  qnrs  J 


(65) 


Given  a  means  for  determining  the  tangential  fields  at  one  of  the  two 
boundaries  and  the  source  distribution,  one  is  able  to  calculate  the 


tangential  fields  at  the  other  boundary  and  therefore  ail  the  fields 
between  the  earth  and  magnetosphere  or  within  the  constant  dipole 
field  cavity  surrounding  the  earth.  Equations  (62)  through  (65)  may  also 


be  written  in  matrix  form, 


EQ(M,q) 

A  i  ...  A  A*  .  ...  A’ 

ql  qN  ql  qN 

Hq(1,  1) 

• 

• 

E.(M,q) 

9 

ql  qN  ql  qN 

H0‘(l,  N) 

H0(Mlq) 

c  ....  c  c* ....  c* 

ql  qN  ql  qN 

H  (1,  1) 

9 

* 

% 

H^(M,q) 

ql  qN  ql  qN 

H  *  (1,  N) 

9  _ 

+  AX  (66) 


a 


iS^SSti 


where 


l  n  •  •  •  Cl  Cl  . 

ql  qN  ql 


lqN  3ql 


•  fl 


qN 


V*’ 


V*0 

Js9(l) 

( 


•  •  §qN  ?'ql  •  '  '  ?qN  ’Tql  '  '  '  'TqN 

•  '  CqN  Cql  '  '  '  CqN  Vql'  '  '  VqN 

— 


I 

Js0<Ni 

P?(D 

PS(N) 


and  q  =  1,  2,  .  .  .  ,  N. 


in.  THE  FORM  OF  THE  FIELDS  AT  THE  OUTER  BOUNDARY 


In  some  manner  the  tangential  fields  at  the  outer  boundary  must 
be  determined  in  order  that  the  matrix  (66)  may  be  inverted  for  the 
solution  of  the  magnetic  field  intensities  at  all  points  on  the  surface  of 
the  inner  spherical  boundary.  The  finite -difference  equations  (38)  through 
(43)  may  then  be  utilized  for  the  final  resolution  of  all  fields  between  the 
two  radial  boundaries. 

Two  alternatives  are  available  for  the  specification  of  the  outer 
boundary  fields.  If  the  scope  of  the  investigation  includes  the  effect  of 
the  boundary  conditions  above  on  the  field  configuration  within  the  con¬ 
fines  of  the  radial  limits,  various  distributions  may  be  assigned  to  the 
tangential  fields  at  the  outer  boundary.  The  problem  then  develops  as 
to  what  physical  significance  the  chosen  boundary  condition  has.  This 
procedure  is  particularly  applicable  to  the  case  in  which  the  exterior 
boundary  is  considered  to  be  a  conductor  supporting  a  surface  current. 

The  alternative  to  this  specification  is  the  mathematical  genera¬ 
tion  of  functions  which  might  conceivably  describe  the  form  of  the  fields 
at  the  outer  boundary.  In  order  that  this  may  be  achieved,  an  investi¬ 
gation  of  the  electromagnetic  waves  produced  in  a  spherical  system  is 


made. 


i 


o 
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The  field  configuration  at  the  outer  limits  of  the  region  of  interest 
is  assumed  to  be  of  the  form  of  spherical  waves  propagating  radilly 
outward.  The  wave  is  considered  to  have  propagated  a  sufficient  distance 
that  the  medium  in  which  the  fields  are  examined  may  be  assumed  to  be 
similar  to  free  space  with  very  small  free  electron  densities. 


A.  Solutions  to  the  Vector  Wave  Equation 

If  the  displacement  currents  are  retained  in  the  Maxwell 


equations,  then 


V  x  H  =  (a  +  iu)e)E  . 


(67) 


The  curl  of  equation  (7)  is 


V  x  V  x  E  =  -  iujVx  (pH) 


(68) 


If  it  is  assumed  that  the  permeability  of  the  medium  is  constant  with 
respect  to  the  space  coordinates,  equation  (68)  becomes 


V  x  V  x  E  =  -iuw-V  x  H , 


(69) 


The  permeability  is  then  taken  to  be  that  of  free  space,  p^ 
Upon  substitution  of  (o7)  into  (69),  we  have 


V  x  V  x  E  =  -  iuup(a  +  iuje)  E » 


or 


where 


VxVxE  =  k  E, 


2  2 

k  =  u)  pe  -  itupo . 


(70) 
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If  the  charge  density  in  the  region  of  interest  is  taken  to  be  negligible, 
equation  (70)  becomes  the  familiar  vector  Helmholtz  equation 

V2E  +  k2E  =  0  (71) 


Except  for  the  simplest  of  cases,  the  vector  Helmholtz 

equation  cannot  be  resolved  into  three  scalar  equations  which  are  easily 

solved.  For  certain  coordinate  systems,  however,  there  exist  three 

•>  — ►  — ► 

independent  vector  solutions  which  are  designated  as  L,  M  and  N  type 
solutions  ^ 


If  n  is  a  constant  vector  with  certain  properties,  then  the 
three  solutions  are 


L  =  V 

(72) 

M  =  V  x  n  ijf 

(73) 

-►  1 

N  =  —  VxVxnty, 

(74) 

where  ^  is  a  function  satisfying  the  scalar  equation 

V2'(i  +  k2  Jr  =  0  (75) 


in  the  particular  coordinate  system  being  used. 

If  we  let  n  =  r  =  r  a^,  the  radial  vector  for  the  spherical 

system,  then 


M 


V  x  ar(r  f) 


f 


M,  when  decomposed  into  its  three  components,  is 


M  =  0 
r 


M  =  -1—  11 
Q  sin  0  30 


m  -  li 
M<f>~  "  ao 
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Likewise, 


(V  x  M)  =  - 


i _  a  .  .  Q  a* . 

r  sin  0  3  0  (SinQ  30 


.  1  lit 

r  sin^  0  dd>2 


<Vx5>e  =  n(H +7*> 

.  _  i  a  a  i 

(  V  x  M) ,  =  — : — r  -r-r  (-r-*  +  —  i|f) 
0  sin  0  30  3r  r  v 


lx  the  vector  operations  of 


VxVxM-k  M  =  0 


are  carried  out,  the  results  are 


(V  x  V  x  M)  -  k  M  =0 
r  r 


(V  x  V  x  M)  ~k2M  _ _ _ _ —  P  — - ‘L(rlli)  + - i - —  (sin  0-^1) 

1  70  M0  sin  0  30L  2  dr[  3r 1  2  .  n  30* W307 

r  r  sm  0 


2  2 

i  a  y  -i  k  _a_ 

2  .  2  2  J  sin  8  30 

r  sin  0  30 


<t> 


(VxVxM)  -k2M  =|5[i?±(r2|±)  +  -Tl— 
r  r  r  r  sin  0 


a©'  a©7 


TTT7  3]  +k2|e('|l)- 

r  sin  0  30 


J 
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But  the  expression  in  each  of  the  brackets  is  just  V^i}i  in  the  spherical 
coordinates.  Therefore, 

(V  x  V  x  M)9  -  k2M0  ^[v2t  +  k2t]  =  0 

(V  x  V  x  M)^  -  k^M  =  —q  +  k%  j  =  0. 

Thus,  since  the  ijf  scalar  is  constructed  to  satisfy  equation  (75),  the 
vector  M  is  a  solution  of  the  vector  equation  (71).  In  a  similar  fashion 
N  may  also  be  shown  to  be  a  solution  of  (71). 

^  —  ► 

The  electric  field  mav  now  be  represented  as  a  M  or  N 
type  field.  It  was  shown  that  M  contains  only  components  transverse 
to  the  radial  direction,  while  N  has  not  only  transverse  components,  but 
also  a  radial  component.  If  the  electric  field  is  described  by  if  =  hi,  3 
becomes 

— ^  1 r  — ► 

H  =  -  A-  Vx  M, 

lUJp 

and  the  system  of  fields  is  transverse  electric  in  nature.  If  E  =  N,  then 
the  fields  are  seen  to  be  transverse  magnetic. 

B.  Derivation  of  the  Function 

For  the  spherical  system  under  consideration,  the  expansion 
of  the  scalar  equation  (75)  is  given  by 


_9_ 

9r 


(r 


2M)+_JL_ 

hr  +  2  . 

r  sin  9 


a  .  .  A  a* 

89(Sm  0  98>  + 


1 


2.2 

r  sin  9 


d<p 


,  2 

k  ijr  = 


(76) 


If  v/e  assume  a  solution  of  the  form 


O  Q 

£-  U 


0*  0)  =  ^r(r)  i'giQ)  ^(0)* 


(77) 


we  have,  after  substitution  of  (77)  into  (76), 


1  d  .  2  d*r 


di 


~  ,  -  - ,  ^  i  i  d  .  .  a  "Ye. 

•if  dr  r  dr  if  sin  6  d0  Sin  d0 

Yr  T  0 


,2 

d  ijj 


<b  .22 
-£  +  k  r  =  0  . 


*<p  sin  0  d0‘ 


(78) 


Because  of  the  separated  state  of  the  variables,  three  equations  are 
extracted  from  equation  (78);  namely. 


2 

d  i; 


d<p 


0  2, 

2  =  -m  ^ 


(79) 


1  d_f  - 

sin  0  d0  Sin  ®  d0  J 


m 


sin  8 
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=  0 


(80) 


.  ci  ik  2 

1  0.2  r  2  p  ,  ,  A 

— r  —  (r  — — )  +  (k  -  ^r)  $  =0, 

2  or  dr  2  Tr 

r  r 


(81) 


2  2 

where  m  and  p  are  separation  constants.  Equation  (80),  when  expanded, 
becomes 


d\  cos0d*8  2  m 

+  :  ~  — ~  +  (p  “ 


d0 


sin  0  d8 


sin  0 


2  . }  ye 


♦  a  =  °- 


After  the  substitution  of  \  =  cos  0,  this  becomes 


dZ  if 


d\ 


1-X 


(82) 
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,  2 

an  equation  having  singularities  at  X.  =  Tl.  If  p  =  n(n  +  l),  the  solution  of 
interest  will  remain  finite  at  these  singularities. 

The  solutions  to  equations  (79)  are  in  the  form  of  sinusoids. 
In  order  that  the  \ji- solutions  be  single- valued,  the  constant  m  must  take 

on  integer  values,  m  =  0,  +1,  +  2, .  The  solution  to  (79)  may  then 

be  written  as 


ijf  =  a  cos  m  cp  +  b  sin  m <p  , 


»  =aeim*  +  be-lm^, 

9 


The  second  form  is  used  in  this  problem  since  this  is  the  variational 
form  assumed  in  the  finite-difference  approach.  If  tne  m  in  equation 
(82)  takes  on  the  value  m  =  0,  there  follows 

2  d  ^'g  dl*Q 

(l-X  ) - Y  -  2X  — [£-  +  n(n  +  l)ty0=  0.  (8 

d\ 

Differentiation  of  (84)  with  respect  to  X  yields 

2  d*  d2^  df  d^ 

(1-x  »  Te<TT>  -  4K  TT-  2  dT+  n<n+1' -3T  *  °- 

aX  dX 


A  second  differentiation  yields 


2  d2  d  ^0  d  d  ^8  d  *n 

(l-X  )— j( — )  -  2(2  +  1)X  — ( - ~-)  +  [n(n+l)  -  2(2+1)] - ~  =  0, 

d\  dX2  dX  dX2  dX2 


and  the  third  time, 
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2  d^  ,  d  d  i|i 

U  ~  ^  )  2  ^ - 3")  "  ~^( - 3~)+  [n(n+l)  -  3(3+1)] - —  =  0. 

dX  dX  dX  dX 
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The  general  expression  for  the  mth  derivative  of  equation  (84)  is 


2  d2  d  *6  d  d 

d),2  d\m  dX  d\m 


,m 

d 

+  [n(n+l)  -  m(m+l)] -  =  0. 

d\m 


Now,  letting 


2 

v  =  (1-X  ) 


equation  (86)  becomes 


?  A2  ?  "m/?  ^  ?~mb 

(1-X  )-^[(l-X  )  v]-2(m+l)X^~  [(1-X  )  v] 

dX 


2  ~mb 

+  [n(n+l)-m(m+l)](l-X  )  v  =  0  .  (87) 


After  a  small  amount  of  algebraic  manipulation,  one  finds  that  (87) reduces 


2  2 

,,  ,2.d  v  ,,  dv  r  ,  m 

(1-X  ) - 5  "  2KdX  +  [n(n+1) - 3  JV  =  °- 

dX  1-X 


Equation  (84)  is  a  special  equation  which  has  solutions  in  the  form  of  the 
Legendre  polynomial  Pn(X),  so  the  solution  to  (88)  must  be 


2m/2  dm 


VC[(1-X)  7^Pn(X>^ 


The  function  within  the  bracket  is  called  the  associated  Legendre  Poly¬ 


nomial  of  order  m  and  degree  n.  It  is  symbolically  noted  as 
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*6=  CPr<X>- 


(89) 


The  third  of  the  three  separated  equations  is 


1  d  2  ^  r  r  2  n(n+l) 

~2  d? <r  -dT)+Ck  °- 

r  r 


The  substitution  of  p  =  kr  and  R  =  tjf^r  leads  to 

1-5:  +  [!- ]  R  =  Oj 

dp  p2 


(90) 


Now,  if  R=  'sTp  S,  equation  (90)  becomes  Bessel*s  equation 


d  S  1  dS 
+ 


dp 


2  p  dp 


+  [i.l£±W2L]s.  o, 

P 


with  solutions  of  the  form  of 


S  =  Zn,l/2<p>> 

where  zn+1/2^p^  represents  the  various  Bessel  functions  of  order  n+1/2. 
The  solution  to  the  separated  equation  in  (81)  is  therefore 

*r=d7-Zn+l/2(p)- 

vr 

The  Hankel  functions  are  chosen  for  this  problem,  and  so 
the  may  be  written  in  terms  of  these  Bessel  functions  of  the  third  kind. 


.  1  „(1)  „  »  1  (2) 

=  d  "  Hn+l/2^kr>  +  6  “  Hn+1/Z<kr>’ 
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and,  using  the  definitions  found  in  Appendix,  the  final  form  of  the  radial 
component  of  the  separated  solutions  is 

=  f  h^  (kr)  +  g  h^  (kr). 

¥r  n  n 

The  combination  of  equations  (83),  (89)  and  (91)  yields  the 
total  form  of  the  solution  to  the  scalar  wave  equation  (75)  in  spherical 
coordinates, 

ill  =  [aeim^  +  be"im<^][cPm(X)][fh(1)(kr)  +  gh(2)(kr)]. 

Two  of  these  terms  may  be  discarded  for  this  problem.  Since  a  depen¬ 
dence  on  <p  of  the  fields  for  the  finite -difference  calculations  was  taken 
to  be  eim^,  the  term  e  irn<^  is  not  required  here. 

The  selection  of  the  proper  Hankel  function  depends  on  the 
fact  that  only  an  outgoing  wave  is  desired.  Now, 

,2  2 

k  =  a)  pe  -  luupa . 

2 

For  a  passive  medium,  a  is  positive  real  so  that  k  will  lie 

in  the  fourth,  and,  for  the  particular  root  chosen,  k  will  have  a  negative 

imaginary  part.  In  order  for  the  wave  to  die  off  at  very  large  distance, 

(2) 

the  second  spherical  Hankel  function,  h^  (kr),  must  be  selected  for 
the  representation  of  the  outward  propagating  wave. 

The  general  solution  is  therefore 

*  t  =  aelm0  h<2)(kr)P™  (\)  (91) 

n  n 

where  the  constant  a  is  a  combination  of  previous  constants. 


C.  The  Final  Form  of  the  Fields 


The  M  and  N  solutions  may  now  be  constructed  using  (91). 


M  =  0 
r 


M 


M 


0 

<t> 


=  jm  eim*h<2)  (kr)pm(x) 
sin  0  “  “ 


n 


n 


dP  (X) 

im <f>  ,  (2)  „  t  n 

=  -  e  h  (kr )  - — - 

n  '  '  d0 


=  £l2±il  elm*h(2)(kr)Pm(k> 

r  kr  n  n 


N  1  eim^  [rh<2)(kr)] 
0  kr  dr  L  n  '  '  J 


dP 


m 


n  (M 


d9 


im 

kr  sin  0 


eim* 


(2) 

n 


(kr)  ]  P 


m 

n 


(M. 


When  \  =  cos  0  is  substituted  into  these  six  equations,  M  and  N  become 


M  =  0 
r 


inn  imd>  ,  (2)„  .  _m,.  . 

M  =  — =-  e  vh'#(kr)P  (X) 

9  «/T\2 

I - dPm(k) 

M .  =  JTTk1  elm*  h  2  (kr) 

"  dX 


<P 


n 


N  =  Otetll  eim,#,h(2,(kr)Pm<k) 
r  kr  n  n  ' 


xt  '/l-X  mn0  d  r  ^(2)  n 

N  =  — - -  e  r-[rh'  (kr)]  — — 

9  kr  dr  n  dX 


dPm(\) 


(92) 

(93) 

(94) 

(95) 

(96) 


N  = - eim<?>  A[rhi  (kr)] pm(X).  (97) 

*  kr*s/I-k2  dr  „ 


Since  any  solenoidal  E  field  may  be  constructed  from  a 
linear  sum  of  transverse  electric  and  transverse  magnetic  waves, 


E  =  V  (A  5  +  BN) ,  (98) 

Z_i  n  n  n  n 

n 

where  and  are  constants  determined  by  the  boundary  conditions. 

In  other  words,  A^  is  the  amplitude  of  the  n**1  mode  transverse  electric 
wave,  and  B^  is  the  amplitude  of  the  nth  rnode  transverse  magnetic  wave. 
The  index  n  must  be  summed  up  to  a  value  high  enough  to  include  the 
highest  ordered  mode  present  in  the  system.  Since  the  finite-difference 
method  employed  N  points  in  the  angular  direction,  the  highest  mode 
we  shall  consider  is  the  N**1  mode. 

The  magnetic  field  H,  as  derived  from  Maxwell's  equation. 


is 


H  =  V  x  {  V  (A  M  +  BN)1, 
nup  k  L>  n  n  n  n  J 


n 


and,  since  A  and  B  are  constant  with  respect  to  the  soace  coordinates, 
n  n 


iouj-t 


I 


n 


K<v 


x  M  )  + 
n 


B  (V  xi  ll  . 
n  n  j 


(99) 


Equation  (74)  defining  N  is  rewritten,  yielding 


N  =  -  V  x  M  , 
k 


(100) 


and  from  equation  (73),  M  is  given  by 
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_►  -*  l  -*  2 

M  =  V  x  Rf  - - j  V  x  n  (-k  \ 

k 


-►  1  -+  2 
M  - - j  V  x  (nV  f)  . 

k^ 


(101) 


From  a  vector  identity  and  the  knowledge  of  the  properties  of  the  vector 
n,  we  have 


V  x  (ViJ  x  n)  =  -  nV  •  Vty  =  -  n  V 


(102) 


Upon  substitution  of  (102)  into  equation  (101),  the  vector 


solution  M  becomes 


M  =  — -  V  x  V  x  (V  x  n) 


It  can  be  seen  that 


M  =  Viji  x  n. 


if  equation  (73)  is  expanded,  and  therefore. 


—  1  „  „  - 
M  =  — r  V  x  V  x  M  , 

k- 


which,  by  equation  (100),  reduces  finally  to 


M  =  -  V  x  N  . 
k 


The  magnetic  field  is  thus  given  most  conveniently  by 


H  =  - 


iiup 


y  [A  N  +  B  M  ] . 
L>  n  n  n  n 


(103) 


n 
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Finally,  with  equations  (92)  through  (97)  substituted  into  (98)  and  (103), 
the  general  field  expressions  for  the  spherical  wave  propagating  out¬ 
wardly  are 

E  (G  2<2+lLeim*h,2,(kr)Pm{\))  (104) 

r  U  \  n  kr  n  n  J 

n=l 


E.  =  y  |f  ^Hl_eim^h‘2)(kr)P>) 

9  L  In  i - -  n  n 


n=l  N/l-X 


-  G  eimj>  -f-  [rh(Z)(kr)]  ■  (105) 

n  kr  dr  n  '  J  dX  J 


dPm(X) 


N  / - -  .  .  dP  ( 

**  ■  l  k^  *imk2w> 


n=l 


+  G  -iH!_  eim*f  [rh(2)(kr)]Pm(M}  (106) 

n  i - t  dr  n  n  J 


kr  Vl-\2 


H  =f  {-r  ^eim<fh<2)(kr)P>)) 

r  Li  l  n  lojur  n  n  J 


n=l 


.m. 


H 


N  •  A  ^  /?X  dPn‘(X) 

r  r„  *1  l-X  im0  d  _  (2)  n 


=  Y  {F  ~ 

Li  l  n  xuul 


,  e  “[rh  (hr)] — — 

8  L  l  n  iu;pr  drL  n  dX 

n=l 


-  G  - =  e 

n 

U)|Jl 


mk  imp  (2) 


n  Jt  \2 
mu.N/1  -X 


h(2)(kr)Pm(X)} 
n  n  J 


(107) 


(108) 


im  <b 
e 


-i[rhJ2)(kr)]P"1(M 

dr  n  n 


"k!  eim^h(2)(kr) 
KJL’p  n 


dPm(X) 

n 

dX 


(109) 


The  four  expressions  of  these  equations  for  the  tangential 
fields  may  be  written  in  symbolic  form  using  notation  introduced  pre- 
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viously. 


N 

E:(M,q)  =  7  [F  a  +  G  a'  ] 
0'  L  n  qn  n  qnJ 


(110) 


n=l 


E  (M,q)  =  )  [F  b  +  G  V  ] 
9  La  n  qn  n  qn 

n=l 


(111) 


H 


N 

,<M,q)=y  [F  f  +G  f*  ] 
3  La  n  qn  n  qn 


n=l 

N 

H  (M,  q)  =7  [F  g  +  G  g‘  ], 
9  La  n°qn  n°qn 

n=l 


(112) 


(113) 


where  the  terms  a  ,  a*  ,  b  ,  b*  ,  f  ,  f1  ,  g  and  g*  are  the 

qn  qn  qn  qn  qn  qn  °qn  qn 

complex  terms  in  the  field  equations  (105),  (106),  (108),  and  (109)  and 
are  functions  of  both  the  mode  number  and  the  angular  position. 


IV.  COMBINATION  OF  THE  EQUATIONS  AND  BOUNDARY  CONDITIONS 


There  now  exist  two  expressions  for  each  of  the  fields  tangential 
to  the  outer  boundary.  Equations  (62)  through  (65)  are  equated  to  (110) 
through  (113),  yielding 


N 

t  [a  F 
L  qn  n 


n=l 


+  a’  G 
qn  n 


+ 


N 

It 

i-i-i 


0L  J 
qn  s<p 


(n)  + 


qn  s <p 


(n>  +  Pqnps(n)]  (H5) 


n=l 


[b  F  +  b*  G  ] 
qn  n  qn  n 


4 


n=l 


[B  H  <1,  n)  +  B'  H  (1.  n)] 
qn  u  qn  <p 


4  tVVn,  +  Vse(n)+VpS(n)]  (U6) 

n=l 


N  N 

V  [f  F  +£'  G  ]  =y  [C  H  (l,n)  +  C'  H  (l.n)] 

lj  qh  n  on  n  on  U  on  (b 


qn  n  qn  n  L>  qi 
n=l  n=l 


4  ^nVn,+  <U7) 

n=l 


N  N 

Y  [g  F  +  g'  G  ]=y  [D  H  (1,  n)  +  D*  H  (1,  n)] 
L  6qn  n  6qn  nJ  L  qn  9  qn  d> 


qn  n  qn  n  L>  qi 
n=l  n=l  ‘ 


N 


+  T  J  *(*)  +  S  J  <n)  +  vnnp  (n)]  (118) 

Li  qn  s <p  qn  so  qn  s 


n=l 


The  quantities  which  are  known  in  these  4N  equations  are  all  the  coeffi¬ 
cients  plus  the  source  terms  J  (n),  J  (n)  and  p  (n)  while  the  unknowns 

SO  So  s 
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are  HA(1,  n),  H  (1,  n),  F  and  G  .  There  are,  therefore,  4  N  equations 
8  <p  n  n 

with  4N  unknowns. 


A.  Evaluation  of  the  Fields  at  the  Inner  Boundary 

Since,  for  each  q,  both  the  coefficients  a  ,  a'  ,  8  ,  v  , 

*  qn  qn  qn  qn 

y',  6  ,  §  ,  V  .  tv  »  C  ,  C'  and  V  and  the  sources  J  (n), 

qn  qn  qn  qn  qn  qn  qn  qn  s<p 

J  (n)  and  p  (n)  are  known  in  the  problem,  let  the  following  substitutions 
s  d  s 

be  made 

Slq  =  -  1  [aqnJs^(n>  +  aqnJs6(n)  +  Sqn'>s(n)  ] 
n=l 


S ,  = 
2q 


S3q  = 


N 

VCY  J  ,  (n)  +  V1  J  Q(n)  +  &  p  (n)  ] 
L  qn  s <f>  ‘qn  s9  qnrs 


n=l 

N 


•  I  [ ?qnJs<J,(n)  +  5qnJ=e(n)  +  n,nPS(n>  ^ 


n=l 


N 


S4q  = 


-  Y  [C  J  An)  +  C'  J  Q(n)  +  V  P  (n)  ] 
L  qn  s<p  bqn  s8  qnrs 


n=l 


Then  equations  (114)  through  (117)  may  be  rewritten  in  matrix  form 

(where  the  coefficients  a„„,  a*  ,  b  ,  f  ,  f*  ,  e  and  g1  are  the 
v  qn’  qn’  qn’  qn  qn  6qn  6qn 

negative  values  of  those  calculated  previously).  This  matrix  is  shown 


on  the  next  page. 


The  system  of  matrices  may  now  be  inverted  by  some  method 


to  give  the  values  of  H.(l,  q)  and  H  (1,  q)  for  each  q,  q  =  1,  2,  .  .  .  ,  N. 

0  < p 

Once  this  is  done,  all  the  tangential  fields  at  the  inner  boundary  are 


determined,  and  the  problem  may  then  be  concluded. 


B.  Field  Calculations 

Once  the  tangential  fields  at  the  inner  boundary  are  completely 

determined,  all  fields  within  the  system  are  calculated  from  the  finite  - 

difference  equations  (38)  through  (43)  in  a  step-by-step  manner  as  before. 

The  proper  boundary  conditions  are  applied  at  the  discontinuity,  and  the 

procedure  is  continued  until  all  six  field  components  are  known  at  each 

of  the  grid  points.  This  concludes  the  procedure  except  for  perhaps  a 

few  checks  on  the  accuracy  of  the  method. 

Although  checks  would  increase  the  length  of  the  computation 

it  might  be  desirable,  at  least  until  the  calculations  prove  to  be  reliable, 

to  perform  one  or  more  checks  on  the  accuracy  of  the  calculations. 

One  such  verification  consists  of  solving  for  the  coefficients 

F  and  G  in  the  matrix  which  was  to  be  inverted.  The  fields  at  the  outer 
n  n 

boundary  calculated  from  those  coefficients  in  equations  (104)  through 
(109)  should  be  the  same  as  those  calculated  in  the  last  iteration  process. 

Any  computational  instability  of  the  fields  either  in  the 
radial  direction  or  polar  angle  direction  must  be  carefully  examined  for 
computational  troubles  such  as  significant  figure  problems. 


V.  COMPUTATIONAL  METHODS 


A  program  for  the  digital  computer  may  be  most  conveniently- 
written  with  several  subprograms  or  subroutines.  A  brief  outline  of 


such  a  program  follow?- 

Quantities  which  must  be  inserted  into  the  computer  initially 

include  the  frequency,  iu;  the  plasma  frequency,  uj  ;  the  number  of 

P 

grid  points  in  the  radial  direction,  M;  the  number  of  grid  points  desired 

for  X,  N;  the  molecular  weight  of  ions  or  mass  ratio  of  ions  to  electrons, 

a  ;  the  radius  of  the  inner  boundary,  R  .  ;  the  radius  of  the  outer 
o  mm 

boundary,  R  ;  the  current  distributions,  JQ(q)  and  J  (q);  the  charge 
max  o  <p 

density,  p(q);  and  the  radial  location  of  the  sources.  Of  course,  con¬ 
stants  such  as  and  tt  must  be  included  in  the  program.  The  quantities 
listed  above  are  read  into  the  computer  to  enable  the  programmer  to 
change  one  or  more  of  these  parameters  easily. 

After  all  inputs  are  within  the  computer,  the  conductivity  matrix 
is  calculated  in  a  separate  subroutine.  This  matrix  has  components  which 
are  found  in  Appendix  B.  These  numbers  are  stored  for  later  use. 

The  next  quantities  to  be  calculated  in  another  subroutine  are  the 

coefficients  of  the  matrix  (66).  The  functions  a  ,  a’  ,  b  ,  b’  ,  f  , 

qn  qn  qn  qn  qn 

P  n»  Sqn»  an<*  §qn  equations  (110)  through  (113)  are  then  computed  with 
separate  subroutines  for  the  calculation  of  the  Hankel  functions  and  the 
associated  Legendre  polynomials.  Once  all  elements  in  the  first  matrix 
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A_  are  determined,  the  matrix  is  inverted  by  Gaussian  elimination  for 

the  field  quantities  HA(1,  q)  and  H,  (1,  q);  this  may  also  be  accomplished 

0  9 

in  a  separate  subroutine. 

The  main  program  may  then  be  used  for  the  final  calculation  of 
all  the  field  quantities.  To  eliminate  the  need  for  even  more  storage 
space,  the  six  field  components  are  printed  as  they  are  calculated  at 


each  radial  level. 


Although  the  use  of  many  subroutines  increases  the  amount  of 
storage  space  required,  it  greatly  simplifies  the  isolation  of  errors  or 
possible  trouble  areas  in  the  program,  as  well  as  simplifying  the  writing 
of  the  program.  If  storage  space  is  the  prime  factor  in  the  calculations, 


as  it  is  likely  to  be,  all  calculations  may  be  accomplished  in  the  main 
program.  If  further  reduction  in  either  the  amount  of  core  storage  or 


time  needed  is  necessary,  the  problem  may  be  solved  in  separate  parts 
or  programs  with  the  output  of  one  program  used  as  the  input  of  the  next. 
In  such  a  way  the  conductivities,  the  coefficients  for  the  matrix  (66),  and 


the  eight  functions  of  (110)  through  (113)  may  be  calculated  separately  and 
used  as  input  to  a  final  program  to  invert  the  matrix  A.,. 

A  high-speed  digital  computer  possessing  an  extremely  large 
memory  capacity  is  needed  for  quantitative  results  from  this  problem. 
Because  of  the  large  storage  volume  needed,  all  calculation  should,  if 
possible,  be  accomplished  using  only  single  precision  arithmetic.  For 
this  reason  means  must  be  found  to  circumvent  any  significant  figure  problems. 
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Several  subtle  difficulties  are  present  in  the  direct  application  of 

some  of  the  difference  equations.  One  such  problem  of  significant  figures 

lies  in  the  calculation  of  from  equation  (30).  For  this  problem  an 

auxiliary  field  E  was  utilized  by  Chapman  (1966)  as  defined  by 
P 


E  =  k  E  +  krtErt 
p  r  r  6  8 


For  further  details  see  reference  3. 


VI.  CONCLUSIONS 


Many  notions  present  themselves  as  to  how  the  quantitative  results 
may  be  used.  One  aspect  of  the  problem  that  certainly  should  be  examined 
is  the  effect  of  various  current  distributions  and  charge  density  distri¬ 
butions  on  the  fields  at  the  surface  of  the  earth.  A  person  need  not  limit 
the  source  to  one  current  sheet,  but,  by  extending  the  idea  of  this  problem, 
he  might  examine  the  fields  with  several  such  current  rings  or  sheets  of 
charge  encircling  the  earth. 

The  variation  of  the  electromagnetic  field  with  the  several  parame¬ 
ters  of  interest  is  another  needed  route  of  investigation.  These  parameters 
include,  among  others,  the  frequency,  plasma  frequency,  altitude  of  the 
prevailing  sources,  and  altitude  of  the  undistorted  dipole  field  cavity. 

In  order  that  an  adequate  resolution  may  be  achieved  and  a  minimum 
chance  of  instability  due  to  computational  errors  be  realized,  a  large 
number  of  grid  points  must  be  involved  in  the  quantitative  consideration 
of  this  problem.  For  this  reason  the  present  digital  computer  available 
(Control  Data  Corporation  1604)  could  only  be  utilized  if  the  problem  were 
divided  into  the  several  parts  mentioned  previously.  Attempts  in  the 
computational  region  will  most  likely  be  made  after  the  installation  of  the 
new  computer  (Control  Data  Corporation  6600)  in  the  Computation  Center 


at  The  University  of  Texas. 


APPENDIX  A 

FINITE-DIFFERENCE  FORMULAE 


The  finite-difference  formulae  are  based  on  polynomial  approxi¬ 
mations.  See  reference  4  for  details  of  the  formulae.  Three  methods 
of  differencing  exist;  namely,  forward  differencing,  backward  differencing, 
and  central  differencing.  For  the  use  of  this  problem,  only  the  expressions 
for  the  derivatives  of  functions  are  of  interest. 

Let  y(x)  be  a  function  of  x  (It  may  also  be  a  function  of  other 
variables,  but  each  of  these  is  held  constant  for  these  approximations 
of  the  partial  derivatives,  and  so  no  generality  is  lost  by  the  consider¬ 
ation  of  only  x),  y^  be  the  value  of  the  function  y  at  the  value  of  x  =  y^f 
y^l  be  the  value  of  the  function  at  the  next  point  chosen  on  the  x-axis 
(x^+^>  x^),  and  y^  ^  be  the  value  of  y  at  the  point  directly  to  the  left  of 

x  on  the  x-axis.  The  figure  illustrates  the  notation, 
n 


Figure  4 

Definitions  for  the  Finite- Difference  Formulae 
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The  forward  difference  formula  for  the  derivative  of  y  with 
respect  to  x  evaluated  at  x  =  is  given  by  the  series 


d y  , 

7n  _  1 

dx  h 


[(yn+ryn>-I<W-2yn+l+yn)  +  -  '  '  ] 


where  h  is  the  increment  on  x  between  x  and  x  The  backward 

n  n+1 

difference  formula  is,  if  h  is  also  the  increment  between  x  ,  and  x  , 

n-1  n 


dy 


dx 


"  =  Kyn‘  Vl>  +  2  (yn‘  2Vl+  yn-2>  '  ‘  ] 


and  the  formula  for  the  more  accurate  central  difference  is 

dy^  ^  j 

dx  2h  L^n+1  ^n-1^  6  ^n+2  2yn+l+2yn-l  yn-2^  +  ’  J* 

If  approximations  to  the  first  order  only  are  used,  the  difference 
formulae  become  simply 


=  yn+l  ~  yn 

dx  x  -  x 
n+1  n 


for  the  forward  difference, 


dy  y  -  y  . 
7n  n  'n-1 


dx  x  -  x 

n  n-1 


for  the  backward  difference,  and 


dyn  _  -  Vl 

^  ’  2(x„+l  -  X  J 


for  the  central  difference.  The  central  difference  formula  is  most  con¬ 


venient  when  the  increments  taken  in  the  x-direction  are  all  equal. 
However,  non-uniform  increments  may  also  be  used,  with  the  various 
formulae  being  modified  accordingly. 

The  central  difference  formula  is  used  whenever  possible,  and 
so  it  is  employed  in  the  computation  of  the  derivatives  of  the  field  v/ith 
respect  to  X  at  all  points.  The  central  difference  is  used  to  approxi¬ 
mate  the  derivatives  with  respect  to  r  except  in  the  initial  step  from 
the  inner  boundary  when  the  forward  difference  is  used. 

Let  AX  be  the  incremental  change  in  X  used  in  the  grid,  and  Ar 
be  the  radial  distance  between  two  adjacent  levels.  Also,  the  notation 
A(p,  q)  is  taken  to  mean  the  value  of  A  at  the  radius  r^  and  angle 
(cosine)X.  A(p+1,  q)  therefore  signifies  the  value  of  A  at  the  next  radial 
and  the  same  value  of  X.  The  difference  formulae  are  then 

dA  _  A(p+1,  q)  -  A(p,  q) 
dr  Ar 

dA  _  A(p,  q)  -  A(p-1,  q) 
dr  Ar 

dA  _  A(p-f3,  q)  -  A(p-1,  q) 
dr  ~  2Ar 

for  the  forward,  backward,  and  central  difference,  respectively,  for 


the  radial  derivative,  and 


i 


\ 


AQ 
A  / 


dA  _  A(p,  q+1)  -  A(p,  q) 
dX.  "  AX. 


dA  _  A(p,  q)  -  A(p,  q-1) 
dX.  "  AX. 


dA  _  A(p,  q+1)  -  A(p,  q-1) 
dX.  "  2  AX. 


for  the  same  differences  of  the  tangential  derivatives.  All  six  of  these 
expressions  are  for  the  derivatives  evaluated  at  the  point  (p,  q). 


APPENDIX  B 


DERIVATION  OF  THE  CONDUCTIVITY  MATRIX 


Bostick  (1964)  derived  the  conductivity  tensor  from  considera¬ 
tions  of  particle  interaction  using  basic  equations  of  plasma  physics. 
The  complex  conductivity  tensor  is  based  on  the  generalized  Ohm's  law, 
the  equation  of  motion  for  the  charged  particles  in  an  ionized  gas,  and 
the  equation  of  motion  for  neutral  particles.  The  following  is  an  outline 
of  the  derivation  of  the  conductivity  elements  used  by  Chapman  (1966). 

It  is  presented  here  since  the  conductivities  must  be  calculated  for  this 
problem. 

Let  f^,  f^  and  f^  be  defined  by 


fl  = 


f,  = 


f,  = 


f. 

f  - 
en  2 

f.  m 

m  e  . 

—5—  +  -  f 

2  m.  en 

l 

m 

f.  +  f  + - f. 

le  en  m.  m 

i 


(A-l) 

(A-2) 

(A- 3) 


where  the  terms  f  ,  f.  and  f.  are  the  effective  collision  frequencies 

en  m  le 

between  electrons  and  neutral  particles,  ions  and  neutral  particles,  and 
ions  and  electrons,  respectively,  and  and  rm  are  the  electron  and 
ion  masses. 

If  the  electron  plasma  frequency  uj  is  defined  as 
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N  q2 
o 

m  e 
e  o 


where  Nq  is  the  electron  charge  density,  then  the  generalized  Ohm's  law 
is  given  by  o 


£(iu>  +  f3)  J  +  U)e  J  x  k  J  =  E  +  Vp  x  kBQ 


e  uj 

o  p 


B  f.  Ap 

+  (V  -  V  )  +  — ^ 

0)  p  n'  N  q 


(A-4) 


In  equation  (A-4),  j  is  the  current  desnity,  k  is  the  unit  vector  in  the  +z 
direction,  ti'e  is  the  cyclotron  resonance  frequency  for  electrons,  is 
the  average  plasma  velocity,  is  the  average  velocity  of  the  neutral 
particles,  and  pg  is  the  perturbation  in  pressure  of  electrons. 

The  equation  of  motion  for  the  charged  component  of  the  plasma 
is 

B  f.  ^  ^ 

iuiNmV  - - —  J  -  N  mf,(V  -V  )  +  jXkB  -  Ap  (A-5) 

o  p  u)  o  2  p  n  o  p 

in  which  m  is  the  ion  mass  and  p^is  the  perturbation  in  pressure 
difference  of  ions  and  electrons. 

The  equation  of  motion  for  the  neutral  particles  is 

B  f  ^ 

iujNmf  = — —  J-N  mfJV  -V  )-Ap  (A- 6) 

n  uj  o  2  n  p'  *n 

e 

where  N  is  the  particle  density  of  neutral  particles,  and  p  is  the 

n 

perturbed  pressure  of  the  neutral  particles. 


Equation  (A-4)  is  rearranged  to  give 


B  f. 

(i ujN  m.  +  N  m.f_)  V  -  N  m.f.  V  =  — - —  J  +  J  x  k B 
oi  oiZpoiZn  u)  o 

e 


(A-7) 


and  equation  (A-6)  is  rearranged  to  give 


-  N  m.f_  V  +  (N  m.f  +  i  ajNm. )  V  = 
oiZp'oiZ  i  n 


B  fi 

o  1 


(A- 8) 


If  equations  (A-7)  and  (A- 8)  are  solved  simultaneously  for  V  and  V  , 

n  p 


the  results  are 


V  = 
P 


TVT  B  £t  B  M, 

. .  .  N  .  .  ol  7  ^  r*_  .  o  1  Z  7 

(f_+i—)( -  J+JxkB  ) - J 

2  N  u)  o  U) 

o  e  e 

- r_ - 

r.  N,  2  N  -I 


N  m.  F i  i»i  f _  (1+  ~)  - 
oil  2  N 


(A- 9) 


V  = 
n 


B  fl  -  B  f, 

— - —  (f  j  +  i  U))  J  +  f ->  ( — - —  J  +  JxkB  ) 
uj  Z  Z  uj  o 

e  e 


(A- 10) 


Then 


V  -V 
P  n 


B  f 

o  1  N  N 

— — —  iiu(l  +  ■—— )  J  +  iour—  B  Jxk 
cu  N  No 

e _ o _ o _ 

r.  ,  .  N  „  Z  N  n 


N  m.  fiu)f,(l+^-)- 
o  i  L  2  N 


(A- 11) 


When  equations  (A-9)  and  (A-ll)  are  substituted  into  equation  (A-4)  and 


the  resulting  terms  rearranged,  there  results 


i  U)  +  fg 

B2 

o 

2 

'  0  U) 
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(A- 12) 


where 


N  m. 


-  (J  x  k)  x  k  [  (i2  +  i  u>  j  =  E 

u  o 


4  . _ T  .2p.  ,  N  .  u)  N  n 
A  =  (N  m.)  [i®f2(l+jj- )  -  T  W  J 

o  2  o 


With  the  cyclotron  frequency  defined  by 

qB. 


m 


then  the  first  term  in  equation  (A-12)  is  reduced  to 

2 

i U)  +  f -  B  f.  iiuN  m. 
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2 - A - (1+iT,J  Vp 
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iiu  +  f^ 
2 

e  u) 
o  p 


f,  m 


V  <— >  r — LJl-^ — 1} « (A-i3) 


e  u>  i  1-Ht-  Vt~,  ) 
op  f  'N  +  N 

^  2  o 


the  second  term  becomes 


Let  a  be  defined  as  the  right-hand  side  of  equation  (A-13),  b  as  the  right- 
hand  side  of  equation  (A-14),  and  c  as  the  right-hand  side  of  equation  (A-15). 
Then 


aj  +  b  J  x  k  +  C(J  xk)xk  =  e  w  E  . 

o  p 


(A-16) 


Now,  the  earth*s  magnetic  field  is  assumed  to  be  of  a  form  pro¬ 
duced  by  a  magnetic  dipole  at  the  center  of  the  earth  with  a  magnetic 
moment  of 


M  =  3.06x10  weber/meters. 


Also,  from  Prince,  Bostick,  and  Smith  (1964),  the  strength  of  the  static 
field  is  given  by 


M  F.  2~“  7T 

=  — =■  n/  4  cos  0  +  sm 


o  3 
r 


or 


_  Mh 
Bo‘  T  ’ 

r 


where 


=  Jl7l\Z  . 


m.  N 

Defining  a  =  —  and  T)  =  -  ,  one  finds  that 

o  m  1  N 
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fj^l+Tl) 

a  =  (£3+1|">'a  ri»+f ,(i+n)] 


b  = 


_  Mqh  f 


m  r 
e 


2f“ 

sl1  a  [i tti  +  f ?(1+Tl)]} 


c  =  -( 


Mqh  ,2  f  +  i  u) 


m  r 
e 


^  {.ittja  [iuu+f  (1+T))]  } 


(A- 17) 


(A-18) 


(A- 19) 


equation  (A-12)  is  specified  in  terms  of  the  rectangular  coordinate 
system,  but  the  problem  for  which  the  conductivity  matrix  is  needed  is  in 
the  frame  o^  the  spherical  coordinate  system.  It  is  therefore  necessary 
that  equation  (A**12)  be  transformed  over  to  the  spherical  system. 


ae  by 


The  unit  vector  k  is  given  in  terms  of  the  unit  vectors  a^  and 


k  =  k  a  +  k  a  , 
r  r  8  0 


where 


k  =  - 
r 


2  cos  0 


n/  1  + 


,  2  A 

3  cos  6 


(A- 20) 


V- 


sin  0 


Jl 


+  3  cos  0 


(A-  21) 


Then  Jxk  becomes 


Jxk  =  (a  J  +  a  J  +a  J  )  x  (k  a  +kaafi) 
rr  80  <j>  <p  rr  0  9 


Jxk  =  -  a  J  k  +  a  J  k  +  a ,  (J  kQ  -  J  k  )  , 
t  <p  &  8  <p  r  0  r9-  0  r 


(A- 22) 


and 


(J  x  k)  x  k  =  a  [k.(J  Jk  -  J  k  J]  +  an  [k  (J  k.  -  Jak  )] 
'  rL  9  8  r  r  9  0 L  r'  r  0  •  0  r,J 


+  V"Vkr  +  k9)] 


2  2 

But  k  +  k_  =  1,  so 
r  0 


(J  x  k)  x  k  =  ar[ke(J6kr  -  Jke)]  +  a8[kr(Jrk0-  J0kr)]  -  .  (A-23) 


Equation  (A-16),  reduced  to  component  form  with  the  substitution  of 


equations  ( A-22)  and  (A-23),  is 


57 


Jr(a-cke2)  +  J0(ckrke)  +  J^l-bkg)  =  e0u)p  (A-24) 

Jr{ckekr)  + J0(a-ckr2)  + J^(bkr)  =  eo.^(  (A-25) 

Jr(bk8)  +  J0(*bkr)  +  J^,(a-C>  =  e0i»pE^  (A-26) 


The  determinant  of  these  three  equations  is 
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=  a[(a  -  c)2+ b2].  (A-27) 
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a  -  c 

The  three  equations  (A-24)  through  (A-26)  are  inverted  to  give  the  current 
distribution  in  terms  of  the  three  electric  field  components.  In  matrix 
form,  the  resulting  equations  are 


<t> 


a  a  Q 

rr  r  0 
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The  resulting  conductivity  matrix  is 
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where  a  is  defined  by  (A-17),  b  by  (A-18),  c  by(A-19),  k  and  k  by 

r  8 

(A- 20)  and  (A-21),  and  D  by  equation  (A-27). 

The  conductivity  matrix  A^  is  very  general,  containing  in  its 
terms  the  effects  of  aprticle  interaction.  If  the  medium  is  assumed  to 
be  such  that  very  few  ion-electron,  ion-neutral  particle,  or  electron- 
neutral  particle  collisions  occur,  elements  contained  in  the  matrix  are 
greatly  simplified. 

For  no  collisions, 

f.  =  f.  =  f  =  0 . 
te  m  en 

Therefore 


and  the  three  terms  a,  b,  and  c  are  found  to  be 


Then  the  various  components  cf  the  matrix  become 


a(a  -  c)  =  -  uu2  [l  — | —  )  1 

ua  u  m  r 
o  e 
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ab  =  iujh 
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f  . 

It  is  seen  that  a(a  -  c)  and  b  -  c(a  -  c)  are  entirely  real  while  ab 


and  D  are  purely  imaginary.  This  causes  a  ,  aQ,,  a  ,  and  a  to 


r 0’  "6<j>  0 r’ 


<p  8 


be  real,  and  a  ,  a  A»  crA  ,  and  a ,  ,  to  be  imaginary.  Also,  it 

r  r  r  o  0r  0  0.  00 

may  be  noted  that  a  =  CJA  ,  a  ,  =  (J,  ,  gn,  =  -ojn. 
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APPENDIX  C 


USEFUL  PROPERTIES  OF  THE  SPHERICAL  BESSEL  FUNCTIONS 
AND  ASSOCIATED  LEGENDRE  POLYNOMIALS 

SPHERICAL  BESSEL  FUNCTIONS 

The  spherical  Bessel  functions  of  the  third  kind  are  defined  by 


h(1)<p)  =  H{1\  (p) 

n  r  2p  n+|  r 


hn  ’<P>  =  ^ 


The  Hankel  functions  H^x  (p)  and  (p)  are  linear  combi- 


n+f 


n+2 


nations  of  the  Bessel  functions  of  the  first  and  second  kins,  namely. 


H(1)  (n\  - 


Jn+l(P>  +  iNn4A<P> 


H<2A(p)  ■  J  ,i(p)  -  iN  i(p) 

n+2  n+2  n+2 


For  large  values  of  p,  the  approximations  for  the  Spherical 
Hankel  functions  are,  as  p -*•«>, 


.(1),  .  1  .  ..n+1  ip 

n  -  p  6 


,  (2)  ,  1  ...n+1  -ip 

n  (p)  -  ^  (l)  e  * 


If  z  (p)  is  either  h^(p)  or  h^(p),  then  some  useful  recurrence 
n  r  n  r  n 
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relations  for  the  spherical  functions  are 


2n  + 1 


P  n 


Z»  =  Zn.i<P)  ‘  Z-xi<P> 


n+1 
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-  [-=-  Z  (p)]  =  (-l)mp‘nZ  (p), 
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for  n  =  0,  +1,  ±  2,  .  .  .  ,  and  m  =  1,  2,  3,  .  .  .  . 


ASSOCIATED  LEGENDRE  POLYNOMIALS 


The  associated  Legendre  polynomial  is  defined  by 


p>>  =  a 
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12 


d\ 


m  n 


P,  (M . 


or 


n  \^\Tn^  jn+m/\  2  un 

pm(X)  =  d  Ql-ZH 

n  ?n  i  ,,n+m 

2  n !  d\ 


Some  recurrence  relations  which  are  useful  in  the  calculation  of 
the  associated  Legendre  polynomials  on  a  digital  computer  are,  with  the 


notation  of  the  functional  dependence  on  X  omitted, 
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